Midterm Assignment

# Dataframe 생성
PM_ratio = pd.DataFrame({"PM10":PM_measured["PM10_indoor"] / PM_measured["PM10_outdoor"], "PM2.5": PM_measured["PM2.5_indoor"] / PM_measured["PM2.5_outdoor"], "PURI": PM_measured["PURI"]})
PM_ratio["dayYear"] = PM_measured["dayYear"]

# 입자 크기별로 나누어 I/O Ratio를 나타낸 그래프
PM_fig2 = px.line(PM_ratio, x = "dayYear", y = ["PM10", "PM2.5"], line_group = "PURI", title = "Historical PM Ratio", labels={"dayYear": "Day", "value": "I/O Ratio (Cin/Cout)", "variable": "Size"})
PM_fig2 = PM_fig2.update_layout(
    font_family="Times New Roman",
    font_color="black",
    font_size = 18,
    title_font_family="Times New Roman",
    title_font_color="black",
    title_font_size = 24,
    legend_title_font_color="black"
)
PM_fig2.write_html("Fig3.html")

# Box plot 확인
fig2 = px.box(PM_ratio, y = ["PM10", "PM2.5"], color = "PURI", labels={"variable": "Size", "value": "I/O Ratio (Cin/Cout)"})
fig2 = fig2.update_layout(
    font_family="Times New Roman",
    font_color="black",
    font_size = 18,
    title_font_family="Times New Roman",
    title_font_color="black",
    title_font_size = 24,
    legend_title_font_color="black"
)
fig2.write_html("Fig4.html")

# 두 데이터의 길이가 다름 -> 샘플링 필요 (무작위로 수행)
PM_ratio_off = PM_ratio[PM_ratio["PURI"] == 0]
PM_ratio_on = PM_ratio[PM_ratio["PURI"] == 1]
PM_ratio_off2 = PM_ratio_off.sample(len(PM_ratio_on), replace = False)

# index 초기화
PM_ratio_off2 = PM_ratio_off2.reset_index()
PM_ratio_on = PM_ratio_on.reset_index()

# 공기 청정기를 껐을 때의 I/O Ratio - 켰을 때의 I/O Ratio
PM_ratio_sub = PM_ratio_off2 - PM_ratio_on
PM_ratio_size = len(PM_ratio_sub)
PM_ratio_df = PM_ratio_size - 1
PM_ratio_mean = PM_ratio_sub.mean()
PM_ratio_se = PM_ratio_sub.std(ddof = 1) / np.sqrt(PM_ratio_size)
print(PM_ratio_mean); print(PM_ratio_se)

# 95% 유의수준에서 신뢰구간 형성
confint_ratio = stats.t.interval(alpha = 0.95, df = PM_ratio_df, loc = PM_ratio_mean[1], scale = PM_ratio_se[1])
confint_ratio2 = stats.t.interval(alpha = 0.95, df = PM_ratio_df, loc = PM_ratio_mean[2], scale = PM_ratio_se[2])

# 두 비의 차이와 0을 비교하는 t-test
PM_ratio_ttest10 = stats.ttest_1samp(PM_ratio_sub["PM10"], 0)
PM_ratio_ttest25 = stats.ttest_1samp(PM_ratio_sub["PM2.5"], 0)
matrix_conf2 = [
    ['Confint Table', 't-value', "Lower", "Upper"],
    ['PM10', PM_ratio_ttest10[0], confint_ratio[0], confint_ratio[1]],
    ['PM2.5', PM_ratio_ttest25[0], confint_ratio2[0], confint_ratio2[1]]
]

# 결과 확인
conf_table2 = ff.create_table(matrix_conf2, index = True)
conf_table2.write_html("Table3.html")

# 데이터 불러오기
data <- read.csv("office_data.csv")
data <- data[complete.cases(data),]

# Cin을 구하는 함수 만들기
pred_Cin <- function(parameters, t, Cin, Cout){
    ACH_P <- parameters[[1]]
    ACH_k <- parameters[[2]]
    S <- parameters[[3]]
    pred <- Cin[1]
    for (i in seq(1, length(t) - 1)){
        h <- t[i+1] - t[i]
        pred[i+1] <- (pred[i] + h / 2 * (ACH_P * (Cout[i] + Cout[i+1]) - ACH_k * pred[i] + 2 * S)) / (1 + h / 2 * ACH_k)
    }
    return (data.frame(dayYear = t, PM_measured = Cin, PM_predicted = pred))
}

# 최소제곱법을 활용하기 위해 모델 예측값과 실측값 사이의 잔차를 구하는 함수 만들기
resid <- function(parameters, t, Cin, Cout){
    temp <- pred_Cin(parameters, t, Cin, Cout)
    return (sum((temp$PM_measured - temp$PM_predicted)^2))
}

# 공기청정기 가동 여부에 따라 데이터를 구분하고, 하루단위로 모드 설정
data_off <- data[which(data$PURI == 0),]
data_on <- data[which(data$PURI == 1),]
data_off$mode <- 0
for (i in seq(14, 1)){
    data_off[which(data_off$dayYear <= 729 + i),]$mode <- i
}
data_on$mode <- 0
for (i in seq(9, 1)){
    data_on[which(data_on$dayYear <= 742 + i),]$mode <- i
}

# Optimization의 Bound 설정
library(DEoptim)
Lp <- c(0, 0, 0)
Up <- c(100, 100, 100)

# PM10, PURI = 0
param1 <- NULL
temp <- NULL
for (i in seq(1, 14)){
    temp <- data_off[which(data_off$mode == i),]
    opt <- DEoptim(resid, Lp, Up, control = DEoptim.control(trace = FALSE, NP = 100, itermax = 500, parallelType = 1, parVar = c("resid", "pred_Cin")), t = temp$dayYear, Cin = temp$PM10_indoor, Cout = temp$PM10_outdoor)
    param1 <- rbind(param1, unlist((opt$optim$bestmem)))
}
param1 <- as.data.frame(param1)
model_data1 <- NULL
for (i in seq(1, 14)){
    model_data1 <- rbind(model_data1, pred_Cin(param1[i,], data_off$dayYear[which(data_off$mode == i)], data_off$PM10_indoor[which(data_off$mode == i)], data_off$PM10_outdoor[which(data_off$mode == i)]))
}

# PM2.5, PURI = 0
param2 <- NULL
temp <- NULL
for (i in seq(1, 14)){
    temp <- data_off[which(data_off$mode == i),]
    opt <- DEoptim(resid, Lp, Up, control = DEoptim.control(trace = FALSE, NP = 100, itermax = 500, parallelType = 1, parVar = c("resid", "pred_Cin")), t = temp$dayYear, Cin = temp$PM2.5_indoor, Cout = temp$PM2.5_outdoor)
    param2 <- rbind(param2, unlist((opt$optim$bestmem)))
}
param2 <- as.data.frame(param2)
model_data2 <- NULL
for (i in seq(1, 14)){
    model_data2 <- rbind(model_data2, pred_Cin(param2[i,], data_off$dayYear[which(data_off$mode == i)], data_off$PM2.5_indoor[which(data_off$mode == i)], data_off$PM2.5_outdoor[which(data_off$mode == i)]))
}

# PM10, PURI = 1
param3 <- NULL
temp <- NULL
for (i in seq(1, 9)){
    temp <- data_on[which(data_on$mode == i),]
    opt <- DEoptim(resid, Lp, Up, control = DEoptim.control(trace = FALSE, NP = 100, itermax = 500, parallelType = 1, parVar = c("resid", "pred_Cin")), t = temp$dayYear, Cin = temp$PM10_indoor, Cout = temp$PM10_outdoor)
    param3 <- rbind(param3, unlist((opt$optim$bestmem)))
}
param3 <- as.data.frame(param3)
model_data3 <- NULL
for (i in seq(1, 9)){
    model_data3 <- rbind(model_data3, pred_Cin(param3[i,], data_on$dayYear[which(data_on$mode == i)], data_on$PM10_indoor[which(data_on$mode == i)], data_on$PM10_outdoor[which(data_on$mode == i)]))
}

# PM2.5, PURI = 1
param4 <- NULL
temp <- NULL
for (i in seq(1, 9)){
    temp <- data_on[which(data_on$mode == i),]
    opt <- DEoptim(resid, Lp, Up, control = DEoptim.control(trace = FALSE, NP = 100, itermax = 500, parallelType = 1, parVar = c("resid", "pred_Cin")), t = temp$dayYear, Cin = temp$PM2.5_indoor, Cout = temp$PM2.5_outdoor)
    param4 <- rbind(param4, unlist((opt$optim$bestmem)))
}
param4 <- as.data.frame(param4)
model_data4 <- NULL
for (i in seq(1, 9)){
    model_data4 <- rbind(model_data4, pred_Cin(param4[i,], data_on$dayYear[which(data_on$mode == i)], data_on$PM2.5_indoor[which(data_on$mode == i)], data_on$PM2.5_outdoor[which(data_on$mode == i)]))
}

param <- rbind(param1, param2, param3, param4)
names(param) <- c("ACH_P", "ACH_k", "S")
param

getRSQ <- function(Exp, Pred){
    rmd=which(is.na(Pred))
    if(length(rmd)>0){
        Exp=Exp[-rmd]
        Pred=Pred[-rmd]
    }
    St <- sum((Exp-mean(Exp))^2)
    S <- sum((Exp-Pred)^2)
    return((St-S)/St)
}
rsq1 <- getRSQ(model_data1$PM_measured, model_data1$PM_predicted) # PM10, off
rsq2 <- getRSQ(model_data2$PM_measured, model_data2$PM_predicted) # PM2.5, off
rsq3 <- getRSQ(model_data3$PM_measured, model_data3$PM_predicted) # PM10, on
rsq4 <- getRSQ(model_data4$PM_measured, model_data4$PM_predicted) # PM2.5, on
rsq <- data.frame("PM10_off" = rsq1, "PM2.5_off" = rsq2, "PM10_on" = rsq3, "PM2.5_on" = rsq4)
rsq

# 모델 예측값들로 데이터프레임 구성
model_data <- rbind(cbind(model_data1, model_data2), cbind(model_data3, model_data4))
names(model_data) <- c("dayYear", "PM10_measured", "PM10_predicted", "dayYear2", "PM2.5_measured", "PM2.5_predicted")
model_data$PURI <- data$PURI

# 그래프 그리기
library(plotly)
fig <- plot_ly(model_data, x = ~dayYear, y = ~PM10_measured, type = "scatter", mode = "lines", name = "PM10_measured", line = list(width = 3, dash = "dot")) %>% 
    add_trace(y = ~PM10_predicted, name = "PM10_predicted", line = list(width = 3, dash = "dash")) %>% 
    add_trace(y = ~PM2.5_predicted, name = "PM2.5_predicted", line = list(width = 3)) %>% 
    add_trace(y = ~PM2.5_measured, name = "PM2.5_measured", line = list(width = 3, dash = "dash")) %>% 
    layout(title = "Comparison of Measured and Predicted Data",
           xaxis = list(title = "Day", zeroline = TRUE),
           yaxis = list(title = "Concentration (ug/m3)"))
fig